#############################################################################################
########### calculate exposure inex for each species ########################################
#############################################################################################

#load packages:
library(R.utils)
library(raster)
library(sp)
library(sf)
library(stringr)
library(circular)
library(plotrix)
library(grDevices)

#Main working dir:
MyDir<-paste("/home/jc217070/Maxent",sep="")  #define working directory on HPC
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_swds")
MUOutDir<-paste(MyDir, "/Outputs2_Crossval/permutation.importance", sep="")
IndSppOutDir<-paste(MyDir, "/Outputs3_Final/Final_Individual_Spp", sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#set project name:
projectName<-"WHOSnakes"

#define 'for' variables:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
timecats<-c("Current","2050_Median","2050_Q90","2050_Q10","2090_Median","2090_Q90","2090_Q10")
#timecats<-c("Current","2050_Median","2090_Median")
#timecats<-c("Current")

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
DatSum<-DatSum[DatSum$Listing_comment!="not listed",]
spp<-DatSum$gen_sp_subtax[DatSum$Listing_comment!="not listed"]
sppA<-DatSum$Accepted_Species[DatSum$Listing_comment!="not listed"]
#spSTART<-which (spp=="Micrurus_mertensi")

 DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
# DatSumMod$RangesizeSQKM_Current[DatSumMod$RangesizeSQKM_Current==0]<-0.01
# DatSumMod$SuitabilitySum_Current[DatSumMod$SuitabilitySum_Current==0]<-0.01
# DatSumMod$Exposure_Current[DatSumMod$Exposure_Current==0]<-0.01
# DatSumMod$Impact_Current[DatSumMod$Impact_Current==0]<-0.01

bg<-raster(paste(PlotDir,"/Cop_FAPAR_mean_2010-2019_WW1km.asc",sep=""), crs="+init=epsg:4326")
myext<-extent(-170,179,-58,90)
bg<-crop(bg, myext)
GLOBraster<-bg
values(GLOBraster)<-ifelse(values(GLOBraster)!=-Inf,0,-Inf) #raster showing all land as 0, ocean as NA

GLOBrasterB<-GLOBraster
values(GLOBrasterB)<-ifelse(values(GLOBrasterB)==0,1,values(GLOBrasterB)) #raster showing all land as 1, ocean as NA
GLOBrasterC<-GLOBraster
GLOBrasterC[GLOBrasterC == 0] <- NA

worldmap  <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))
SDMregions   <- st_read(paste(PlotDir,"/worldmap/SDMRegions_WGS84.shp",sep=""))

#make region data frame columns
########################################################################################################################################################################################################################
#basic info:
myregions <- unique(SDMregions$Short_Name)     #vector of region names
myregionsLong <- unique(SDMregions$Region_Nam)     #vector of region names
# region.ras <- data.frame(myregions,myregionsLong) #empty vector to add data for regions into
# names(region.ras)<-c("region_name_short","region_name_long")
# region.ras$n_spp                       <- NA
# ########################################################################################################################################################################################################################
# #columns for range size changes:
# #2050
# region.ras$n_RangesizeSQKM_contraction5to50P_2050_Median<-NA 
# region.ras$n_RangesizeSQKM_expansion5to50P_2050_Median <-NA
# region.ras$n_RangesizeSQKM_contraction50plusP_2050_Median <-NA
# region.ras$n_RangesizeSQKM_expansion50plusP_2050_Median<-NA
# region.ras$n_RangesizeSQKM_stable_2050_Median<-NA
# #2090
# region.ras$n_RangesizeSQKM_contraction5to50P_2090_Median<-NA
# region.ras$n_RangesizeSQKM_expansion5to50P_2090_Median <-NA
# region.ras$n_RangesizeSQKM_contraction50plusP_2090_Median <-NA
# region.ras$n_RangesizeSQKM_expansion50plusP_2090_Median <-NA
# region.ras$n_RangesizeSQKM_stable_2090_Median <-NA
# ########################################################################################################################################################################################################################
# #columns for suitability changes:
# #2050
# region.ras$n_SuitabilitySum_contraction5to50P_2050_Median<-NA
# region.ras$n_SuitabilitySum_expansion5to50P_2050_Median<-NA
# region.ras$n_SuitabilitySum_contraction50plusP_2050_Median<-NA
# region.ras$n_SuitabilitySum_expansion50plusP_2050_Median <-NA
# region.ras$n_SuitabilitySum_stable_2050_Median <-NA
# #2090
# region.ras$n_SuitabilitySum_contraction5to50P_2090_Median<-NA
# region.ras$n_SuitabilitySum_expansion5to50P_2090_Median<-NA
# region.ras$n_SuitabilitySum_contraction50plusP_2090_Median <-NA
# region.ras$n_SuitabilitySum_expansion50plusP_2090_Median<-NA
# region.ras$n_SuitabilitySum_stable_2090_Median<-NA
# ########################################################################################################################################################################################################################
# #columns for exposure changes:
# #2050
# region.ras$n_Exposure_contraction5to50P_2050_Median<-NA
# region.ras$n_Exposure_expansion5to50P_2050_Median<-NA
# region.ras$n_Exposure_contraction50plusP_2050_Median<-NA
# region.ras$n_Exposure_expansion50plusP_2050_Median<-NA
# region.ras$n_Exposure_stable_2050_Median<-NA
# #2090
# region.ras$n_Exposure_contraction5to50P_2090_Median <-NA 
# region.ras$n_Exposure_expansion5to50P_2090_Median<-NA
# region.ras$n_Exposure_contraction50plusP_2090_Median <-NA
# region.ras$n_Exposure_expansion50plusP_2090_Median<-NA
# region.ras$n_Exposure_stable_2090_Median<-NA
# ########################################################################################################################################################################################################################
# #columns for shift vectors:
# #2050
# region.ras$med_ShiftDist_2050_Median<-NA
# region.ras$med_ShiftDir_2050_Median<-NA
# region.ras$mean_ShiftDist_2050_Median<-NA
# region.ras$mean_ShiftDir_2050_Median<-NA
# region.ras$q10_ShiftDist_2050_Median<-NA
# region.ras$q10_ShiftDir_2050_Median <-NA
# region.ras$q90_ShiftDist_2050_Median<-NA
# region.ras$q90_ShiftDir_2050_Median <-NA
# region.ras$med_ShiftDist_2050_Q10 <-NA
# region.ras$med_ShiftDir_2050_Q10 <-NA
# region.ras$med_ShiftDist_2050_Q90 <-NA
# region.ras$med_ShiftDir_2050_Q90<-NA
# #2090
# region.ras$med_ShiftDist_2090_Median<-NA
# region.ras$med_ShiftDir_2090_Median<-NA
# region.ras$mean_ShiftDist_2090_Median<-NA
# region.ras$mean_ShiftDir_2090_Median<-NA
# region.ras$q10_ShiftDist_2090_Median<-NA
# region.ras$q10_ShiftDir_2090_Median <-NA
# region.ras$q90_ShiftDist_2090_Median<-NA
# region.ras$q90_ShiftDir_2090_Median <-NA
# region.ras$med_ShiftDist_2090_Q10 <-NA
# region.ras$med_ShiftDir_2090_Q10<-NA 
# region.ras$med_ShiftDist_2090_Q90 <-NA 
# region.ras$med_ShiftDir_2090_Q90<-NA
# ########################################################################################################################################################################################################################
# 
# 
# 
# 
# for (region in myregions){
#  
#   print(paste("filling out data table for region",region))
#   n<- which(region.ras$region_name_short==region)
#   regionpoly <- subset(SDMregions, SDMregions$Short_Name == region)     #polygon for this region
#   sf_use_s2(FALSE)
#   mypoliesTF<-st_intersects(WHOpolies,regionpoly,sparse=FALSE) #species polygons overlapping with this region polygon
#   WHOpolies$TF<-mypoliesTF
#   mypolies<-subset(WHOpolies, WHOpolies$TF=="TRUE")
#   myspp<- sort(gsub("_"," ",unique(mypolies$Scientific)))    #vector of species present in this region
#   myspp<- myspp[myspp %in% gsub("_"," ",DatSum$gen_sp_subtax)] #safety, don't use any polygons for species that are not listed
#   region.ras$n_spp[n] <- length(myspp)
#   
#   ########################################################################################################################################################################################################################
#   #columns for range size changes:
#   #2050
#   region.ras$n_RangesizeSQKM_contraction5to50P_2050_Median[n]    <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   <  (-0.05) &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current >  (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_RangesizeSQKM_expansion5to50P_2050_Median[n]      <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   >    0.05  &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current <    0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_RangesizeSQKM_contraction50plusP_2050_Median[n]   <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   <= (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_RangesizeSQKM_expansion50plusP_2050_Median[n]     <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   >=   0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_RangesizeSQKM_stable_2050_Median[n]               <- length(DatSumMod$RangesizeSQKM_CHANGE_2050_Median[DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current   <=    0.05  &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2050_Median/DatSumMod$RangesizeSQKM_Current >=  (-0.05) &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   #2090
#   region.ras$n_RangesizeSQKM_contraction5to50P_2090_Median[n]    <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   <  (-0.05) &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current >  (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_RangesizeSQKM_expansion5to50P_2090_Median[n]      <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   >    0.05  &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current <    0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_RangesizeSQKM_contraction50plusP_2090_Median[n]   <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   <= (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_RangesizeSQKM_expansion50plusP_2090_Median[n]     <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   >=   0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_RangesizeSQKM_stable_2090_Median[n]               <- length(DatSumMod$RangesizeSQKM_CHANGE_2090_Median[DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current   <=    0.05  &
#                                                                                                                         DatSumMod$RangesizeSQKM_CHANGE_2090_Median/DatSumMod$RangesizeSQKM_Current >=  (-0.05) &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   ########################################################################################################################################################################################################################
#   #columns for suitability changes:
#   #2050
#   region.ras$n_SuitabilitySum_contraction5to50P_2050_Median[n]    <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   <  (-0.05) &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   >  (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_SuitabilitySum_expansion5to50P_2050_Median[n]      <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   >    0.05  &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   <    0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_SuitabilitySum_contraction50plusP_2050_Median[n]   <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   <= (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_SuitabilitySum_expansion50plusP_2050_Median[n]     <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   >=   0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_SuitabilitySum_stable_2050_Median[n]               <- length(DatSumMod$SuitabilitySum_CHANGE_2050_Median[DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   <=    0.05  &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2050_Median/DatSumMod$SuitabilitySum_Current   >=  (-0.05) &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   #2090
#   region.ras$n_SuitabilitySum_contraction5to50P_2090_Median[n]    <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   <  (-0.05) &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   >  (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_SuitabilitySum_expansion5to50P_2090_Median[n]      <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   >    0.05  &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   <    0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_SuitabilitySum_contraction50plusP_2090_Median[n]   <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   <= (-0.5)  &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_SuitabilitySum_expansion50plusP_2090_Median[n]     <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   >=   0.5   &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_SuitabilitySum_stable_2090_Median[n]               <- length(DatSumMod$SuitabilitySum_CHANGE_2090_Median[DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   <=    0.05  &
#                                                                                                                         DatSumMod$SuitabilitySum_CHANGE_2090_Median/DatSumMod$SuitabilitySum_Current   >=  (-0.05) &
#                                                                                                                         DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   ########################################################################################################################################################################################################################
#   #columns for exposure changes:
#   #2050
#   region.ras$n_Exposure_contraction5to50P_2050_Median[n]    <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        <  (-0.05) &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current          >  (-0.5)  &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_Exposure_expansion5to50P_2050_Median[n]      <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        >    0.05  &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current          <    0.5   &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_Exposure_contraction50plusP_2050_Median[n]   <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        <= (-0.5)  &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_Exposure_expansion50plusP_2050_Median[n]     <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        >=   0.5   &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_Exposure_stable_2050_Median[n]               <- length(DatSumMod$Exposure_CHANGE_2050_Median[DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current                        <=    0.05  &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2050_Median/DatSumMod$Exposure_Current          >=  (-0.05) &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   #2090
#   region.ras$n_Exposure_contraction5to50P_2090_Median[n]    <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        <  (-0.05) &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current          >  (-0.5)  &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_Exposure_expansion5to50P_2090_Median[n]      <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        >    0.05  &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current          <    0.5   &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_Exposure_contraction50plusP_2090_Median[n]   <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        <= (-0.5)  &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_Exposure_expansion50plusP_2090_Median[n]     <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        >=   0.5   &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   region.ras$n_Exposure_stable_2090_Median[n]               <- length(DatSumMod$Exposure_CHANGE_2090_Median[DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current                        <=    0.05  &
#                                                                                                                           DatSumMod$Exposure_CHANGE_2090_Median/DatSumMod$Exposure_Current          >=  (-0.05) &
#                                                                                                                           DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)])
#   ########################################################################################################################################################################################################################
#   #columns for shift vectors:
#   #2050
#   region.ras$med_ShiftDist_2050_Median[n]            <- median(DatSumMod$DistM_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.ras$med_ShiftDir_2050_Median[n]             <- median.circular(circular(DatSumMod$DirDeg_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.ras$mean_ShiftDist_2050_Median[n]            <- mean(DatSumMod$DistM_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.ras$mean_ShiftDir_2050_Median[n]             <- mean.circular(circular(DatSumMod$DirDeg_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                 units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.ras$q10_ShiftDist_2050_Median[n]            <- quantile(DatSumMod$DistM_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],probs=0.1, na.rm=TRUE)
#   region.ras$q10_ShiftDir_2050_Median[n]             <- quantile.circular(circular(DatSumMod$DirDeg_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                               units = "degrees",modulo = "asis",rotation = "clock"),probs=0.1, na.rm=TRUE)
#   region.ras$q90_ShiftDist_2050_Median[n]            <- quantile(DatSumMod$DistM_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],probs=0.9, na.rm=TRUE)
#   region.ras$q90_ShiftDir_2050_Median[n]             <- quantile.circular(circular(DatSumMod$DirDeg_2050_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                               units = "degrees",modulo = "asis",rotation = "clock"), probs=0.9,na.rm=TRUE)
#   region.ras$med_ShiftDist_2050_Q10[n]               <- median(DatSumMod$DistM_2050_Q10[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.ras$med_ShiftDir_2050_Q10[n]                <- median.circular(circular(DatSumMod$DirDeg_2050_Q10[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.ras$med_ShiftDist_2050_Q90[n]               <- median(DatSumMod$DistM_2050_Q90[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.ras$med_ShiftDir_2050_Q90[n]                <- median.circular(circular(DatSumMod$DirDeg_2050_Q90[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   #2090
#   region.ras$med_ShiftDist_2090_Median[n]            <- median(DatSumMod$DistM_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.ras$med_ShiftDir_2090_Median[n]             <- median.circular(circular(DatSumMod$DirDeg_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.ras$mean_ShiftDist_2090_Median[n]            <- mean(DatSumMod$DistM_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.ras$mean_ShiftDir_2090_Median[n]             <- mean.circular(circular(DatSumMod$DirDeg_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.ras$q10_ShiftDist_2090_Median[n]            <- quantile(DatSumMod$DistM_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],probs=0.1, na.rm=TRUE)
#   region.ras$q10_ShiftDir_2090_Median[n]             <- quantile.circular(circular(DatSumMod$DirDeg_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                    units = "degrees",modulo = "asis",rotation = "clock"), probs=0.1,na.rm=TRUE)
#   region.ras$q90_ShiftDist_2090_Median[n]            <- quantile(DatSumMod$DistM_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],probs=0.9, na.rm=TRUE)
#   region.ras$q90_ShiftDir_2090_Median[n]             <- quantile.circular(circular(DatSumMod$DirDeg_2090_Median[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                    units = "degrees",modulo = "asis",rotation = "clock"), probs=0.9,na.rm=TRUE)
#   region.ras$med_ShiftDist_2090_Q10[n]               <- median(DatSumMod$DistM_2090_Q10[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.ras$med_ShiftDir_2090_Q10[n]                <- median.circular(circular(DatSumMod$DirDeg_2090_Q10[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   region.ras$med_ShiftDist_2090_Q90[n]               <- median(DatSumMod$DistM_2090_Q90[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)], na.rm=TRUE)
#   region.ras$med_ShiftDir_2090_Q90[n]                <- median.circular(circular(DatSumMod$DirDeg_2090_Q90[DatSumMod$Accepted_Species %in% gsub("_"," ",myspp)],
#                                                                                  units = "degrees",modulo = "asis",rotation = "clock"), na.rm=TRUE)
#   ########################################################################################################################################################################################################################
# }
# 
# write.table(x=region.ras, file=paste(MyDir, "/Maxent_LUTables/",projectName,"_RegionData.csv",sep=""),
#             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
# 
# 
# 
# ######################
# #plot bar graphs for regions:
# 
 #rearrange order of regions for plot:
 myregions<-c(myregions[3],myregions[4],myregions[2],myregions[10],myregions[6],
              myregions[7],myregions[1],myregions[8],myregions[5],myregions[9])
# 
# png(filename=paste(PlotDir,"/Region_CCBars.png",sep=""), units="in", width=16, height=23, res=500) #png smaller file size, A3 for better resolution
# par(mfrow = c(5,2), mar = c(5,5,5,1), oma=c(0,0.5,2,0.5), mgp=c(1.7,0.5,0))
# # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
# 
# 
# for (region in myregions){
#   
#   n<- which(region.ras$region_name_short==region)
#   
#   # Create matrix
#   set.seed(112)
#   data <- matrix(nrow=5,ncol=6)
#   colnames(data) <- c("Range Size 2050","Range Size 2090","Suitability 2050","Suitability 2090","Exposure 2050","Exposure 2090")
#   rownames(data) <- c("<-50%","<-5%","+/-5%",">+5%",">+50%")
#   
#   # Get data for this region:
#   data[,"Range Size 2050"]<-c(region.ras$n_RangesizeSQKM_contraction50plusP_2050_Median[n],
#                               region.ras$n_RangesizeSQKM_contraction5to50P_2050_Median[n],
#                               region.ras$n_RangesizeSQKM_stable_2050_Median[n],
#                               region.ras$n_RangesizeSQKM_expansion5to50P_2050_Median[n],
#                               region.ras$n_RangesizeSQKM_expansion50plusP_2050_Median[n])/region.ras$n_spp[n]
#   data[,"Range Size 2090"]<-c(region.ras$n_RangesizeSQKM_contraction50plusP_2090_Median[n],
#                               region.ras$n_RangesizeSQKM_contraction5to50P_2090_Median[n],
#                               region.ras$n_RangesizeSQKM_stable_2090_Median[n],
#                               region.ras$n_RangesizeSQKM_expansion5to50P_2090_Median[n],
#                               region.ras$n_RangesizeSQKM_expansion50plusP_2090_Median[n])/region.ras$n_spp[n]
#   data[,"Suitability 2050"]<-c(region.ras$n_SuitabilitySum_contraction50plusP_2050_Median[n],
#                               region.ras$n_SuitabilitySum_contraction5to50P_2050_Median[n],
#                               region.ras$n_SuitabilitySum_stable_2050_Median[n],
#                               region.ras$n_SuitabilitySum_expansion5to50P_2050_Median[n],
#                               region.ras$n_SuitabilitySum_expansion50plusP_2050_Median[n])/region.ras$n_spp[n]
#   data[,"Suitability 2090"]<-c(region.ras$n_SuitabilitySum_contraction50plusP_2090_Median[n],
#                               region.ras$n_SuitabilitySum_contraction5to50P_2090_Median[n],
#                               region.ras$n_SuitabilitySum_stable_2090_Median[n],
#                               region.ras$n_SuitabilitySum_expansion5to50P_2090_Median[n],
#                               region.ras$n_SuitabilitySum_expansion50plusP_2090_Median[n])/region.ras$n_spp[n]
#   data[,"Exposure 2050"]<-c(region.ras$n_Exposure_contraction50plusP_2050_Median[n],
#                                region.ras$n_Exposure_contraction5to50P_2050_Median[n],
#                                region.ras$n_Exposure_stable_2050_Median[n],
#                                region.ras$n_Exposure_expansion5to50P_2050_Median[n],
#                                region.ras$n_Exposure_expansion50plusP_2050_Median[n])/region.ras$n_spp[n]
#   data[,"Exposure 2090"]<-c(region.ras$n_Exposure_contraction50plusP_2090_Median[n],
#                                region.ras$n_Exposure_contraction5to50P_2090_Median[n],
#                                region.ras$n_Exposure_stable_2090_Median[n],
#                                region.ras$n_Exposure_expansion5to50P_2090_Median[n],
#                                region.ras$n_Exposure_expansion50plusP_2090_Median[n])/region.ras$n_spp[n]
#   
#   
#   # Get the stacked barplot
#   barplot(data, 
#           col=c("royalblue","lightblue1","grey90","pink","red"), 
#           border="white", 
#           space=c(0.1,0.1,0.4,0.1,0.4,0.1), 
#           ylim=c(0,1),
#           font.axis=1, cex.axis=1.2,
#           xlab=NULL,ylab="Number of Species",cex.lab=1.2,
#           names.arg=c("2050","2090","2050","2090","2050","2090"),
#           legend.tex=FALSE,
#           main=region.ras$region_name_long[n],cex.main=1.8)
#   #args.legend = list(x = "topright",inset = c(-0.2, 0),cex=0.8)
#   lines (c(0,7.3),c(0.5,0.5), col="grey", lwd= 2, lty=1)
#   
#   #add text for range size, suitability, etc below x axis:
#   # text(0.5,0.5, labels=paste("Range Size"),pos=4, col="black", cex=0.5)
#   # text(3,-1.5, labels=paste("Suitability"),pos=4, col="black", cex=0.5)
#   # text(5.5,-1.5, labels=paste("Exposure"),pos=4, col="black", cex=0.5)
#   
#   mtext(text=paste("Range Size"), font= 1, side = 1, line = 2.2, cex= 1,adj=0.135)
#   mtext(text=paste("Suitability"), font= 1, side = 1, line = 2.2, cex= 1,adj=0.5)
#   mtext(text=paste("Exposure"), font= 1, side = 1, line = 2.2, cex= 1,adj=0.87)
#   
#   # text(0.05, data[1,1]/2, labels=paste(">50%\ndecrease"),pos=4, col="black", cex=0.9)
#   # text(0.05, data[1,1]+data[2,1]/2, labels=paste("5% to 50%\ndecrease"),pos=4, col="black", cex=0.9)
#   # text(0.05, data[1,1]+data[2,1]+data[3,1]/2, labels=paste("stable\n(0% to 5%\nchange)"),pos=4, col="black", cex=0.9)
#   # text(0.05, data[1,1]+data[2,1]+data[3,1]+data[4,1]/2, labels=paste("5% to 50%\nincrease"),pos=4, col="black", cex=0.9)
#   # text(0.05, data[1,1]+data[2,1]+data[3,1]+data[4,1]+data[5,1]/2, labels=paste(">50%\nincrease"),pos=4, col="black", cex=0.9)
# 
# 
# }
# 
# dev.off()
# 
# #############################
# #plot circular graphs for regions
# 
# 
# 
# #############################
# #plot map with species vectors, and region vectors
# 
# mycols1<-c("grey90","grey85","grey80","grey75","grey70",
#            "grey65","grey60","grey55","grey50","grey45",
#            "grey40","grey35","grey30")
# xs <- c(24,55,46,133,78,-58,114,-101,102,-99)
# ys <- c(8,40,58,-25,33,-12,0,21,39,45)
# 
# 
# png(filename=paste(PlotDir,"/regionshiftmap.png",sep=""), units="in", width=16, height=16, res=500) #png smaller file size, A3 for better resolution
# par(mfrow = c(2,1), mar = c(1,0,1,0), oma=c(1,1,1,0), mgp=c(1.2,0.5,0))
# #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
# 
# #Plot 1: 2050
# #plot background
# plot(bg, col=mycols1, box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE,main=paste("2050"))
# plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
# plot(SDMregions[2], border="black", lwd= 2, add=TRUE,col=NA,reset=FALSE)
# 
# #plot species vectors
# for (sp in spp){
#   l<-which (DatSumMod$Accepted_Species==gsub("_"," ",spp))
#   arrows(x0=DatSumMod$long_Current[l], 
#          y0=DatSumMod$lat_Current[l], 
#          x1=DatSumMod$long_2050_Median[l], 
#          y1=DatSumMod$lat_2050_Median[l],
#          length= 0.08, angle=20,col="yellow3",lty=1,lwd=1.2)
# }
# 
# for (region in myregions){
#   
#   n<- which(region.ras$region_name_short==region)
#   s<- which(myregions==region)
#   
#   DirDeg<-    c(region.ras$med_ShiftDir_2050_Q10[n],
#                 region.ras$med_ShiftDir_2050_Q90[n],
#                 region.ras$med_ShiftDir_2050_Median[n],
#                 region.ras$mean_ShiftDir_2050_Median[n])
#   DirDist<-   c(region.ras$med_ShiftDist_2050_Q10[n],
#                 region.ras$med_ShiftDist_2050_Q90[n],
#                 region.ras$med_ShiftDist_2050_Median[n],
#                 region.ras$mean_ShiftDist_2050_Median[n])/1000/100*20 # meters div by 1000=km, div by 100=dd
#   DirDegCirc<- circular(DirDeg, type = "directions", 
#                         units = "degrees",
#                         modulo = "asis", 
#                         rotation = "clock")
#   #Draw arrows with package circular:
#   arrows.circular(DirDegCirc,DirDist, x0 = rep(xs[s],4), y0 = rep(ys[s],4), na.rm = FALSE,zero = (pi/2),
#                   col=c("royalblue","orangered","black","black"),
#                   lwd=c(2,2,2,2),lty=c(1,1,1,3))
# }
# 
# 
# 
# #Plot 2: 2090
# #plot background
# plot(bg, col=mycols1,box=FALSE, bty="n", cex.axis=0.8, legend=FALSE, axes=FALSE,main=paste("2090"))
# plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
# plot(SDMregions[2], border="black", lwd= 2, add=TRUE,col=NA,reset=FALSE)
# 
# #plot species vectors
# for (sp in spp){
#   l<-which (DatSumMod$Accepted_Species==gsub("_"," ",spp))
#   arrows(x0=DatSumMod$long_Current[l], 
#          y0=DatSumMod$lat_Current[l], 
#          x1=DatSumMod$long_2090_Median[l], 
#          y1=DatSumMod$lat_2090_Median[l],
#          length= 0.08, angle=20,col="yellow3",lty=1,lwd=1.2)
# }
# 
# for (region in myregions){
#   
#   n<- which(region.ras$region_name_short==region)
#   s<- which(myregions==region)
#   
#   DirDeg<-    c(region.ras$med_ShiftDir_2090_Q10[n],
#                 region.ras$med_ShiftDir_2090_Q90[n],
#                 region.ras$med_ShiftDir_2090_Median[n],
#                 region.ras$mean_ShiftDir_2090_Median[n])
#   DirDist<-   c(region.ras$med_ShiftDist_2090_Q10[n],
#                 region.ras$med_ShiftDist_2090_Q90[n],
#                 region.ras$med_ShiftDist_2090_Median[n],
#                 region.ras$mean_ShiftDist_2090_Median[n])/1000/100*20 # meters div by 1000=km, div by 100=dd
#   DirDegCirc<- circular(DirDeg, type = "directions", 
#                       units = "degrees",
#                       modulo = "asis", 
#                       rotation = "clock")
#   #Draw arrows with package circular:
#     arrows.circular(DirDegCirc,DirDist, x0 = rep(xs[s],4), y0 = rep(ys[s],4), na.rm = FALSE,zero = (pi/2),
#                     col=c("royalblue","orangered","black","black"),
#                     lwd=c(2,2,2,2),lty=c(1,1,1,3))
# }
# 
# dev.off()

##############################################################
# show main 10 species per region that have the highest current exposure index
# 
 DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
 region.ras<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_RegionData.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
# 
# png(filename=paste(PlotDir,"/Region_Main10Spp.png",sep=""), units="in", width=16, height=23, res=500) #png smaller file size, A3 for better resolution
# par(mfrow = c(5,2), mar = c(10,10,5,1), oma=c(0,0.5,2,0.5), mgp=c(5,0.5,0))
# #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
# 
# for (region in myregions){
#   
#   n<- which(region.ras$region_name_short==region)
#   s<- which(myregions==region)
#   
#   print(paste(Sys.time(), "extracting exposure data per species for", myregions[s], sep=" "))
#   
#   regionpoly <- subset(SDMregions, SDMregions$Short_Name == region)     #polygon for this region
#   regionXs<-st_coordinates(regionpoly)[,1]
#   regionYs<-st_coordinates(regionpoly)[,2]
#   myextMOD<-extent(regionpoly)
#   snapraster<-crop(GLOBrasterC, myextMOD)
#   regionras  <- rasterize(x=regionpoly,y=snapraster,as.numeric(1))
#   
#   sf_use_s2(FALSE)
#   mypoliesTF<-st_intersects(WHOpolies,regionpoly,sparse=FALSE) #species polygons overlapping with this region polygon
#   WHOpolies$TF<-mypoliesTF
#   mypolies<-subset(WHOpolies, WHOpolies$TF=="TRUE")
#   myspp<- sort(gsub("_"," ",unique(mypolies$Scientific)))    #vector of species present in this region
#   myspp<- myspp[myspp %in% gsub("_"," ",DatSum$gen_sp_subtax)] #safety, don't use any polygons for species that are not listed
# 
#   sp.exps<-vector()
#   
#   for (sp in myspp){
#     
#     print(paste("calculating exposure for",sp,"within",region))
#     
#     outfile<-paste(IndSppOutDir,"/",gsub(" ","_",sp),"_Current_Exp.asc",sep="")
#     
#     # unzip raster
#     if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
#       gunzip(filename=paste(outfile,".gz", sep=""),
#              destname=paste(outfile),
#              remove=FALSE, overwrite=TRUE)
#     }
#     
#     # read raster:
#     Exp<-raster(paste(outfile), crs="+init=epsg:4326")
#     #only keep raster values within region:
#     Exp<-Exp*regionras
#     
#     myvals <- values(Exp)
#     myvals<- myvals[myvals>0]
#     myvals<- myvals[!is.na(myvals)]
#     mysum<-sum(myvals)
# 
#     sp.exps<-c(sp.exps,mysum)
#   
#     #remove unzipped raster or zip
#     if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
#       file.remove(outfile)
#     }
#     if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
#       gzip(filename=paste(outfile))
#     }
#   
#   }
# 
#   # DatSumModRegion<-DatSumMod[DatSumMod$Accepted_Species %in% myspp,]
#   DatSumModRegion<-data.frame(myspp,sp.exps)
#   names(DatSumModRegion)<-c("Accepted_Species","Exposure_Current")
#   DatSumModRegion<-DatSumModRegion[order(DatSumModRegion$Exposure_Current,decreasing=TRUE),]
# 
#   print(paste(Sys.time(),"exposure for",region,":"))
#   print(paste(DatSumModRegion))
#   
#   if(s==1){
#     RegionExps<-DatSumModRegion
#   } else{
#     RegionExps<-rbind(RegionExps,DatSumModRegion)
#   }
#   
#   sp.exps<-DatSumModRegion$Exposure_Current[1:10]
#   tot.exp<-sum(DatSumModRegion$Exposure_Current,na.rm=TRUE)
#   res.exp<-tot.exp-sum(sp.exps,na.rm=TRUE)
#   sp.exps<-c(sp.exps,res.exp)
#   top10<-DatSumModRegion$Accepted_Species[1:10]
#   top10<-gsub(" ","\n",c(top10,'other'))
#   
#   barplot(sp.exps,las=2,
#           col=c(heat.colors(11,1,rev=FALSE)), 
#           border="white", 
#           space=c(0.1), 
#           ylim=c(0,max(DatSumMod$Exposure_Current)),
#           font.axis=1, cex.axis=1.2,
#           xlab=NULL,ylab="EI",cex.lab=1.2,
#           names.arg=top10,
#           legend.tex=FALSE,
#           main=region.ras$region_name_long[n],cex.main=1.8)
# 
#   # polygon(regionpoly, border = NA,
#   #         col=rgb(col2rgb("black")[1],col2rgb("black")[2],col2rgb("black")[3],max=255,alpha=75))
#   
# }
# 
# dev.off()
# 
# write.table(x=RegionExps, file=paste(MyDir, "/Maxent_LUTables/PerBioRegionExposures.csv",sep=""),
#             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
# 
# 
############################
#do for each country

mycountries <- gsub("C<U+00F4>te D<U+2019>Ivoire","Cote d'Ivoire",worldmap$CODE_CTY)   #vector of region names

for (country in mycountries){

  n<- which(mycountries==country)
  countryName<- gsub("C<U+00F4>te D<U+2019>Ivoire","Cote d'Ivoire",worldmap$NAME_CTY[n])
  
  print(paste(Sys.time(), "extracting exposure data per species for", countryName, sep=" "))

  regionpoly <- subset(worldmap, worldmap$CODE_CTY == country)     #polygon for this region
  # regionXs<-st_coordinates(regionpoly)[,1]
  # regionYs<-st_coordinates(regionpoly)[,2]
  myextMOD<-extent(regionpoly)
  
  sf_use_s2(FALSE)
  mypoliesTF<-st_intersects(WHOpolies,regionpoly,sparse=FALSE) #species polygons overlapping with this region polygon
  WHOpolies$TF<-mypoliesTF
  mypolies<-subset(WHOpolies, WHOpolies$TF=="TRUE")
  myspp<- sort(gsub("_"," ",unique(mypolies$Scientific)))    #vector of species present in this region
  myspp<- myspp[myspp %in% gsub("_"," ",sppA)] #safety, don't use any polygons for species that are not listed

  if(length(myspp)>=1){
    
  snapraster<-crop(GLOBrasterC, myextMOD)
  regionras  <- rasterize(x=regionpoly,y=snapraster,as.numeric(1))
    
  sp.exps<-vector()

  for (sp in myspp){

    print(paste(Sys.time(),"calculating exposure for",sp,"within",country))

    outfile<-paste(IndSppOutDir,"/",gsub(" ","_",sp),"_Current_Exp.asc",sep="")

    # unzip raster
    if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
      gunzip(filename=paste(outfile,".gz", sep=""),
             destname=paste(outfile),
             remove=FALSE, overwrite=TRUE)
    }

    # read raster:
    Exp<-raster(paste(outfile), crs="+init=epsg:4326")
    #only keep raster values within region:
    Exp<-Exp*regionras

    myvals <- values(Exp)
    myvals<- myvals[myvals>0]
    myvals<- myvals[!is.na(myvals)]
    mysum<-sum(myvals)

    sp.exps<-c(sp.exps,mysum)

    #remove unzipped raster or zip
    if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
      file.remove(outfile)
    }
    if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
      gzip(filename=paste(outfile))
    }

  }

  # DatSumModRegion<-DatSumMod[DatSumMod$Accepted_Species %in% myspp,]
  DatSumModRegion<-data.frame(myspp,sp.exps)
  names(DatSumModRegion)<-c("Accepted_Species","Exposure_Current")
  DatSumModRegion<-DatSumModRegion[order(DatSumModRegion$Exposure_Current,decreasing=TRUE),]
  DatSumModRegion$CountryCode<-paste(country)
  DatSumModRegion$CountryName<-paste(countryName)
  
  print(paste(Sys.time(),"exposure for",country,":"))
  print(paste(DatSumModRegion))

  if(n==1){
    CountryExps<-DatSumModRegion
  } else{
    CountryExps<-read.table(file=paste(MyDir, "/Maxent_LUTables/PerCountryExposures.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
    CountryExps<-rbind(CountryExps,DatSumModRegion)
  }

  write.table(x=CountryExps, file=paste(MyDir, "/Maxent_LUTables/PerCountryExposures.csv",sep=""),
              na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
  
  sp.exps<-DatSumModRegion$Exposure_Current[1:10]
  tot.exp<-sum(DatSumModRegion$Exposure_Current,na.rm=TRUE)
  res.exp<-tot.exp-sum(sp.exps,na.rm=TRUE)
  sp.exps<-c(sp.exps,res.exp)
  top10<-DatSumModRegion$Accepted_Species[1:10]
  top10<-gsub(" ","\n",c(top10,'other'))

  png(filename=paste(PlotDir,"/",country,"Main10Spp.png",sep=""), units="in", width=16, height=16, res=500) #png smaller file size, A3 for better resolution
  par(mfrow = c(1,1), mar = c(10,10,5,1), oma=c(0,0.5,2,0.5), mgp=c(5,0.5,0))
  # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
  
  barplot(sp.exps,las=2,
          col=c(heat.colors(11,1,rev=FALSE)),
          border="white",
          space=c(0.1),
          ylim=c(0,max(DatSumModRegion$Exposure_Current)),
          font.axis=1, cex.axis=1.2,
          xlab=NULL,ylab="EI",cex.lab=1.2,
          names.arg=top10,
          legend.tex=FALSE,
          main=countryName,cex.main=1.8)

  # polygon(regionpoly, border = NA,
  #         col=rgb(col2rgb("black")[1],col2rgb("black")[2],col2rgb("black")[3],max=255,alpha=75))

  dev.off()
  
  removeTmpFiles(h=0.1) #removes all temporary files older than 30min, then read hotspot map back in to continue process
  
  }else{
    
    print(paste(Sys.time(),"no snakes present in",country,":"))
    
  }

}

write.table(x=CountryExps, file=paste(MyDir, "/Maxent_LUTables/PerCountryExposures.csv",sep=""),
            na = "NA", row.names = FALSE, col.names = TRUE, sep=",")


######################################################

# DatSumMod<-read.table(file=paste(MyDir, "/Maxent_LUTables/",projectName,"_SupMatFinFin.csv",sep=""), header=TRUE, sep=",") #load data summary for all species
# 
# DatSumMod$noEDRorSDM<-NA
# DatSumMod$EDRonly<-NA
# DatSumMod$SDMonly<-NA
# DatSumMod$EDRandSDM<-NA
# DatSumMod$SDMsuitOUT<-NA
# DatSumMod$SDMsuitIN<-NA
# DatSumMod$SDMsuitOUTmax<-NA
# 
# 
# #for each species {
# for (sp in spp){
#   
#   print(paste(Sys.time(),"calculating EDR & SDM overlap for",sp))
#   
#   n<-which(DatSumMod$Accepted_Species==gsub("_"," ",sp))
#   
#   #read in SDM
#   outfile<-paste(IndSppOutDir,"/",sp,"_Current_CropThreshCDCut.asc",sep="")
#   
#   # unzip raster
#     if(file.exists(paste(outfile,".gz", sep="")) & !file.exists(outfile)){
#     gunzip(filename=paste(outfile,".gz", sep=""),
#            destname=paste(outfile),
#            remove=FALSE, overwrite=TRUE)
#     }
# 
#     # read raster
#     SDM.ras<-raster(paste(outfile), crs="+init=epsg:4326")
# 
#     #read in polygon
#     mypolies <- subset(WHOpolies, WHOpolies$Scientific == gsub("_"," ",sp))
#     #rasterize polygon (10)
#     snapraster<-crop(GLOBrasterC, extent(SDM))
#     mypolies<-as_Spatial(mypolies)
#     EDR.ras  <- rasterize(x=mypolies,y=snapraster,as.numeric(10))
#     #make NA=0
#     EDR.ras[is.na(EDR.ras)] <- 0
# 
#     #add SDM to polygon 
#     my.ras<-SDM.ras+EDR.ras
#     
#     #remove unzipped file
#     if(file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
#             file.remove(outfile)
#           }
#           if(!file.exists(paste(outfile,".gz", sep="")) & file.exists(outfile)){
#             gzip(filename=paste(outfile))
#           }
# 
#         # get all values
#     myvals <- values(my.ras)
#     myvals<- myvals[!is.na(myvals)]
#     
#     #get all values & calculate length of each (x=0: both say absent, 0<x<1: not in EDR but in SDM, 
#     #                                           x=10: only in EDR but NOT SDM, 10<x<=11: both in EDR & SDM)
#     myvalsNoNo<- myvals[myvals==0]
#     myvalsEDRyes<- myvals[myvals==10]
#     myvalsSDMyes<- myvals[myvals>0 & myvals<=1]
#     myvalsYesYes<- myvals[myvals>10]
#     
#     DatSumMod$noEDRorSDM[n]<-length(myvalsNoNo)
#     DatSumMod$EDRonly[n]<-length(myvalsEDRyes)
#     DatSumMod$SDMonly[n]<-length(myvalsSDMyes)
#     DatSumMod$EDRandSDM[n]<-length(myvalsYesYes)
# 
#     #calculate sum of all 0<x<1 (total suitability predicted outside EDR)
#     #and max (shows if any high suitability predictions lie outside EDR)
#     #and subtract 10 from all 10<x<=11 & calculate sum (total suitability predicted inside EDR)
#     DatSumMod$SDMsuitOUT[n]<-sum(myvalsSDMyes)
#     DatSumMod$SDMsuitOUTmax[n]<-max(myvalsSDMyes,na.rm=TRUE)
#     myvalsYesYes<-myvalsYesYes-10
#     DatSumMod$SDMsuitIN[n]<-sum(myvalsYesYes)
#     
# }
# 
# DatSumMod$SDMtotal<- DatSumMod$SDMonly+DatSumMod$EDRandSDM
# DatSumMod$SDMsuittotal<-DatSumMod$SDMsuitOUT+DatSumMod$SDMsuitIN
# 
# #save summary table
# write.table(x=DatSumMod, file=paste(MyDir, "/Maxent_LUTables/PerCountryExposures.csv",sep=""),
#             na = "NA", row.names = FALSE, col.names = TRUE, sep=",")
# 
# 
# png(filename=paste(PlotDir,"/EDRSDMboxplots.png",sep=""), units="in", width=8, height=16, res=500) #png smaller file size, A3 for better resolution
# par(mfrow = c(2,1), mar = c(10,10,5,1), oma=c(0,0.5,2,0.5), mgp=c(5,0.5,0))
# # #mar=c(bottom, left, top, right); oma= outer margins; mpg = axis label locations c(3, 1, 0)
# 
# #1. plot boxplot of areas (length) inside and outside EDR
# boxplot(list(
#   DatSumMod$SDMonly/SDMtotal,
#   DatSumMod$EDRandSDM/SDMtotal,
#   DatSumMod$EDRonly/SDMtotal,
#   DatSumMod$SDMsuitOUT/SDMsuittotal,
#   DatSumMod$SDMsuitIN/SDMsuittotal,
#   ylab=paste("Fraction of total SDM value (0-1)"),
#   names=c("Commission\nArea","Congruence\nArea","Omission\nArea","Congruence\nSuitability", "Commission\nSuitability"),
#   col=c("white","grey90","grey50","grey30","black")
# ))
# 
# #commission=false positive, omission=false negative
# 
# #2. plot boxplot of total suitability (length) inside and outside EDR, and max outside EDR
# boxplot(list(
#   DatSumMod$SDMsuitOUT/SDMtotal,
#   DatSumMod$SDMsuitIN/SDMtotal,
#   DatSumMod$SDMsuitOUTmax,
#   ylab=paste("Suitability (0-1)"),
#   names=c("Average\nCongruence\nSuitability","Average\nCommission\nSuitability","Maximum\nCommission\nSuitability"),
#   col=c("grey90","grey40","black")
# ))


#############################################################################################
#############################################################################################
#############################################################################################